Historical Renewables Era#

This tutorial explores the Historical_2023_etys scenario using the full ETYS network topology (2000+ buses). By 2023, Great Britain’s electricity system had undergone significant transformation with coal phase-out complete and wind/solar providing substantial generation.

What You’ll Learn#

  • Working with the detailed ETYS transmission network

  • Analyzing high renewable penetration scenarios

  • Identifying transmission constraints and congestion hotspots

  • Understanding renewable curtailment patterns

  • Evaluating system flexibility requirements

Prerequisites#

Run the workflow to generate the solved network:

snakemake resources/network/Historical_2023_etys_solved.nc -j 4

Scenario Overview#

Parameter

Value

Modelled Year

2023

Network Model

ETYS (2000+ buses)

Data Source

DUKES (thermal), REPD (renewables), ESPENI (demand)

Solve Period

First week of January

Key Features

No coal, high wind/solar, network detail

[1]:
# Import required libraries
import os
import sys
import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
from plotly.subplots import make_subplots
from pyproj import Transformer
from pathlib import Path

# Configure matplotlib and plotly
plt.rcParams.update({'font.size': 12})
plt.style.use('ggplot')
%matplotlib inline

# Configure Plotly to output HTML for Sphinx/nbsphinx compatibility
pio.renderers.default = 'notebook_connected'

# Suppress warnings
import warnings
warnings.filterwarnings('ignore')

print('✓ Libraries imported successfully')
✓ Libraries imported successfully
[2]:
# Load solved network (resolve path robustly from different working directories)
network_rel = Path('resources/network/Historical_2023_etys_solved.nc')
network_path = network_rel
# If the relative path doesn't exist from the current cwd, search upward for the repository root
if not network_path.exists():
    for parent in [Path.cwd()] + list(Path.cwd().parents)[:6]:
        candidate = parent / network_rel
        if candidate.exists():
            network_path = candidate
            break
    else:
        # leave as relative; we'll raise a clear error below if not found
        network_path = network_rel

print(f'Loading network from: {network_path}')

if not network_path.exists():
    raise FileNotFoundError(f'Network file not found. Tried: {network_path}.\\nMake sure "resources/network/Historical_2023_etys_solved.nc" exists relative to the repository root.')

network = pypsa.Network(str(network_path))

# Derive a scenario id from the network filename when not provided
scenario_id = 'Historical_2023_etys'

print(f'✓ Network loaded successfully')
print(f'  Buses: {len(network.buses)}')
print(f'  Generators: {len(network.generators)}')
print(f'  Lines: {len(network.lines)}')
print(f'  Storage units: {len(network.storage_units)}')
print(f'  Timesteps: {len(network.snapshots)}')
print(f'  Period: {network.snapshots[0]} to {network.snapshots[-1]}')
Loading network from: c:\Users\alyden\OneDrive - University of Edinburgh\Python\PyPSA-GB v0.0.1\resources\network\Historical_2023_etys_solved.nc
INFO:pypsa.network.io:Imported network 'Historical_2023_etys (Full)' has buses, carriers, generators, lines, links, loads, storage_units, sub_networks, transformers
✓ Network loaded successfully
  Buses: 2044
  Generators: 4766
  Lines: 1592
  Storage units: 81
  Timesteps: 168
  Period: 2023-05-15 00:00:00 to 2023-05-21 23:00:00

Network Summary#

[3]:
# Calculate generation by carrier
p_by_carrier = network.generators_t.p.groupby(
    network.generators.carrier, axis=1).sum()

# Add storage output
storage_by_carrier = network.storage_units_t.p.groupby(
    network.storage_units.carrier, axis=1).sum()
storage_by_carrier[storage_by_carrier < 0] = 0

p_by_carrier = pd.concat([p_by_carrier, storage_by_carrier], axis=1)

# Add H2_turbine links (hydrogen power plants modeled as links, not generators)
if len(network.links) > 0:
    h2_turbine_links = network.links[network.links['carrier'] == 'H2_turbine']
    if len(h2_turbine_links) > 0 and len(network.links_t.p1) > 0:
        # p1 is the power output at bus1 (GB side), which is negative
        h2_output = network.links_t.p1[h2_turbine_links.index].abs()
        h2_ts = pd.DataFrame({'H2_turbine': h2_output.sum(axis=1)})
        p_by_carrier = pd.concat([p_by_carrier, h2_ts], axis=1)
        print(f'Added H2_turbine generation: {h2_turbine_links.p_nom.sum():.0f} MW capacity')

# Separate interconnectors from internal power lines
# Interconnectors connect to external buses (e.g., 'HVDC_External_France')
if len(network.links) > 0:
    # Identify interconnectors vs internal links by checking for 'external' in bus names
    is_interconnector = (
        network.links.bus0.str.lower().str.contains('external', na=False) |
        network.links.bus1.str.lower().str.contains('external', na=False)
    )
    interconnector_links = network.links[is_interconnector].index
    internal_links = network.links[~is_interconnector].index

    # Get interconnector imports (positive p0 = flow into GB from external)
    if len(interconnector_links) > 0:
        ic_flows = network.links_t.p0[interconnector_links].copy()
        ic_flows[ic_flows < 0] = 0
        interconnector_import = pd.DataFrame({'Interconnectors Import': ic_flows.sum(axis=1)})
        p_by_carrier = pd.concat([p_by_carrier, interconnector_import], axis=1)

    # Report internal links separately (not included in generation mix)
    if len(internal_links) > 0:
        print(f'Note: {len(internal_links)} internal power transmission links (not shown in generation mix)')
else:
    print('No links in network')

print('Generation Mix (MWh):')
print(p_by_carrier.sum().sort_values(ascending=False))
print()
print(f'Total demand satisfied: {network.loads_t.p_set.sum().sum():,.0f} MWh')
Note: 3 internal power transmission links (not shown in generation mix)
Generation Mix (MWh):
wind_offshore                      1.150376e+06
EU_import                          9.750260e+05
wind_onshore                       6.977331e+05
nuclear                            5.686804e+05
CCGT                               4.091246e+05
waste_to_energy                    3.105649e+05
solar_pv                           2.760712e+05
coal                               1.719459e+05
Interconnectors Import             1.631320e+05
landfill_gas                       1.024287e+05
Pumped Storage Hydroelectricity    4.560986e+04
biogas                             4.497526e+04
advanced_biofuel                   4.392902e+04
OCGT                               3.257013e+04
large_hydro                        2.311478e+04
Battery                            2.195361e+04
sewage_gas                         7.551726e+03
small_hydro                        5.620384e+03
shoreline_wave                     2.765148e+03
tidal_stream                       6.709302e+01
load_shedding                      6.912225e+00
dtype: float64

Total demand satisfied: 4,805,482 MWh

Interactive Generation Mix Visualization#

[4]:
# Interactive stacked area chart with Plotly
# Define distinct colors for different technologies (ordered: baseload -> renewables -> peakers)
color_map = {
    # Baseload (bottom of stack)
    'nuclear': '#9467BD',          # Purple - distinctive for nuclear
    # Low-carbon dispatchable
    'waste_to_energy': '#8C564B',  # Brown
    'landfill_gas': '#D2691E',     # Chocolate
    'biogas': '#2E8B57',           # Sea green
    'sewage_gas': '#556B2F',       # Dark olive
    'advanced_biofuel': '#6B8E23', # Olive drab
    'biomass': '#228B22',          # Forest green
    'Other Bioenergy': '#3CB371',  # Medium sea green (grouped)
    # Renewables (variable)
    'wind_offshore': '#1E90FF',    # Dodger blue
    'wind_onshore': '#00CED1',     # Dark turquoise
    'solar_pv': '#FFD700',         # Gold
    'large_hydro': '#4169E1',      # Royal blue
    'small_hydro': '#87CEEB',      # Sky blue
    'tidal_stream': '#20B2AA',     # Light sea green
    'shoreline_wave': '#48D1CC',   # Medium turquoise
    'Other Renewables': '#5F9EA0', # Cadet blue (grouped)
    # Storage (mid-merit)
    'battery': '#32CD32',          # Lime green
    'Pumped Storage Hydroelectricity': '#00FA9A', # Medium spring green
    'pumped_hydro': '#00FA9A',
    # Gas (flexible)
    'CCGT': '#FF6347',             # Tomato red
    'OCGT': '#FF4500',             # Orange red (peaker)
    # Coal and other fossil
    'coal': '#696969',             # Dim grey
    # Imports
    'Interconnectors Import': '#DDA0DD', # Plum
    'EU_import': '#DA70D6',        # Orchid
    # Hydrogen system
    'H2': '#00FFFF',               # Cyan (legacy H2 generators)
    'H2_turbine': '#00CED1',       # Dark turquoise (H2 to power)
    'electrolysis': '#40E0D0',     # Turquoise (power to H2)
    'H2_gas': '#7FFFD4',           # Aquamarine (H2 storage)
    # Emergency
    'load_shedding': '#DC143C'     # Crimson
}

# Define stacking order: baseload at bottom, peakers at top
stacking_order = [
    'nuclear',                    # Baseload (bottom)
    'waste_to_energy', 'landfill_gas', 'biogas', 'sewage_gas', 'advanced_biofuel', 'biomass', 'Other Bioenergy',
    'wind_offshore', 'wind_onshore', 'solar_pv',  # Renewables
    'large_hydro', 'small_hydro', 'tidal_stream', 'shoreline_wave', 'Other Renewables',
    'battery', 'Pumped Storage Hydroelectricity', 'pumped_hydro',  # Storage
    'H2', 'H2_turbine',           # Hydrogen power
    'CCGT',                       # Mid-merit gas
    'Interconnectors Import', 'EU_import',  # Imports
    'coal',                       # Coal
    'OCGT',                       # Peaker
    'load_shedding'               # Emergency (top)
]

# Group small contributors (<1% of total) into 'Other' categories
total_gen = p_by_carrier.sum().sum()
threshold = total_gen * 0.01  # 1% threshold

# Identify small bioenergy and renewable types to group
bioenergy_types = ['biogas', 'sewage_gas', 'advanced_biofuel', 'landfill_gas']
renewable_types = ['tidal_stream', 'shoreline_wave', 'small_hydro']

p_by_carrier_grouped = p_by_carrier.copy()

# Group small bioenergy
small_bio = [c for c in bioenergy_types if c in p_by_carrier.columns and p_by_carrier[c].sum() < threshold]
if len(small_bio) > 0:
    p_by_carrier_grouped['Other Bioenergy'] = p_by_carrier_grouped[small_bio].sum(axis=1)
    p_by_carrier_grouped = p_by_carrier_grouped.drop(columns=small_bio)

# Group small renewables
small_ren = [c for c in renewable_types if c in p_by_carrier.columns and p_by_carrier[c].sum() < threshold]
if len(small_ren) > 0:
    p_by_carrier_grouped['Other Renewables'] = p_by_carrier_grouped[small_ren].sum(axis=1)
    p_by_carrier_grouped = p_by_carrier_grouped.drop(columns=small_ren)

# Select columns with generation > 0
cols_with_gen = p_by_carrier_grouped.sum()[p_by_carrier_grouped.sum() > 0].index.tolist()

# Sort by stacking order
cols_to_plot = [c for c in stacking_order if c in cols_with_gen]
# Add any remaining columns not in stacking_order
cols_to_plot += [c for c in cols_with_gen if c not in cols_to_plot]

# Create interactive stacked area chart
fig = go.Figure()

for col in cols_to_plot:  # Stack in order
    fig.add_trace(go.Scatter(
        x=p_by_carrier_grouped.index,
        y=p_by_carrier_grouped[col] / 1e3,  # Convert to GW
        name=col.replace('_', ' ').title(),
        mode='lines',
        stackgroup='one',
        fillcolor=color_map.get(col, '#CCCCCC'),
        line=dict(width=0.5, color=color_map.get(col, '#CCCCCC')),
        hovertemplate='<b>%{fullData.name}</b><br>' +
                      'Time: %{x}<br>' +
                      'Power: %{y:.2f} GW<br>' +
                      '<extra></extra>'
    ))

fig.update_layout(
    title=dict(text='Generation Mix - Historical_2023_etys', x=0.5, xanchor='center'),
    xaxis_title='Time',
    yaxis_title='Generation (GW)',
    hovermode='x unified',
    height=600,
    legend=dict(
        orientation='v',
        yanchor='top',
        y=1,
        xanchor='left',
        x=1.02
    ),
    template='plotly_white'
)

fig.show()
print('Stacking order: Baseload (bottom) -> Renewables -> Storage -> Gas -> Peakers (top)')
print('Tip: Click legend items to show/hide technologies. Double-click to isolate one.')
Stacking order: Baseload (bottom) -> Renewables -> Storage -> Gas -> Peakers (top)
Tip: Click legend items to show/hide technologies. Double-click to isolate one.

Interactive Network Topology#

Explore the network topology with interactive maps showing bus locations, transmission lines/links, and generators overlaid on a GB map.

Understanding Network Topology#

The maps below show the physical layout of Great Britain’s electricity infrastructure:

  • Buses: Connection points (substations, generation sites) shown as orange markers

  • Lines/Links: Transmission corridors carrying power between buses

  • Geography: Real GB coordinates using British National Grid projection

The network model captures how electricity flows are constrained by transmission capacity, creating locational price differences and potential congestion.

[5]:
# Interactive network map using PyPSA's explore() function
# Select a representative snapshot for visualization
snapshot_idx = len(network.snapshots) // 2
snapshot = network.snapshots[snapshot_idx]

print(f'Visualizing network at snapshot: {snapshot}')

# Detect network type based on bus count
is_zonal = len(network.buses) < 50 and len(network.lines) == 0 and len(network.links) > 0
is_reduced = len(network.buses) < 100 and len(network.buses) >= 30
network_type = 'Zonal' if is_zonal else ('Reduced' if is_reduced else 'ETYS')
print(f'Network type detected: {network_type} ({len(network.buses)} buses)')

# Network topology visualization (electricity only, no hydrogen)
try:
    import copy
    import folium
    if 'x' in getattr(network.buses, 'columns', []) and 'y' in getattr(network.buses, 'columns', []):
        # Filter out hydrogen buses and links
        h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas', 'H2 pipeline']
        h2_buses = set()
        if 'carrier' in network.buses.columns:
            h2_buses = set(network.buses[network.buses.carrier.isin(h2_carriers)].index)

        # Create filtered network copy
        network_plot = copy.deepcopy(network)

        # Remove hydrogen buses if any found
        if len(h2_buses) > 0:
            for bus_id in h2_buses:
                if bus_id in network_plot.buses.index:
                    network_plot.remove('Bus', bus_id)

        # Also filter hydrogen links
        if len(network_plot.links) > 0 and 'carrier' in network_plot.links.columns:
            h2_links = network_plot.links[network_plot.links.carrier.isin(h2_carriers)].index
            for link_id in h2_links:
                if link_id in network_plot.links.index:
                    network_plot.remove('Link', link_id)

        # Convert OSGB -> WGS84 lon/lat for mapping
        t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
        lon, lat = t.transform(network_plot.buses['x'].to_numpy(), network_plot.buses['y'].to_numpy())
        network_plot.buses['x'] = lon
        network_plot.buses['y'] = lat

        print(f'Electricity network: {len(network_plot.buses)} buses, {len(network_plot.lines)} lines, {len(network_plot.links)} links')

        # Create interactive map
        m = network_plot.plot.explore(map_style='light', tooltip=True)
        display(m)
        print('✓ Network topology displayed (electricity only)')
    else:
        print('⚠️  No bus coordinates (x, y) found. Network visualization requires coordinate data.')
except Exception as e:
    print(f'⚠️  Network topology map unavailable: {e}')
    import traceback
    traceback.print_exc()
Visualizing network at snapshot: 2023-05-18 12:00:00
Network type detected: ETYS (2044 buses)
Electricity network: 2044 buses, 1592 lines, 13 links
✓ Network topology displayed (electricity only)
[6]:
# Interactive network map showing transmission loading using PyPSA explore
try:
    import copy

    import folium
    from folium import plugins

    # Detect transmission element type (electricity only)
    has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
    has_internal_links = False

    if has_lines:
        # ETYS/Reduced network - use lines
        loading = (network.lines_t.p0.loc[snapshot].abs() / network.lines.s_nom).fillna(0)
        component_type = 'lines'
    elif len(network.links) > 0:
        # Zonal network - use internal electricity links only
        h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas', 'H2 pipeline']
        is_external = (
            network.links.bus0.str.lower().str.contains('external', na=False) |
            network.links.bus1.str.lower().str.contains('external', na=False)
        )
        is_h2_link = network.links.carrier.isin(h2_carriers)
        internal_links = network.links[~is_external & ~is_h2_link]
        if len(internal_links) > 0 and snapshot in network.links_t.p0.index:
            link_flows = network.links_t.p0.loc[snapshot, internal_links.index].abs()
            loading = (link_flows / internal_links.p_nom).fillna(0)
            component_type = 'links'
            has_internal_links = True
        else:
            loading = pd.Series(dtype=float)

    if len(loading) > 0:
        # Convert bus coordinates
        t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
        bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
        bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)

        # Create folium map
        center_lat = bus_coords['lat'].mean()
        center_lon = bus_coords['lon'].mean()
        m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')

        # Color function
        def get_color_by_loading(load_val):
            if load_val >= 0.95: return '#DC143C'  # Crimson (>95%)
            elif load_val >= 0.80: return '#FF8C00'  # Orange (80-95%)
            elif load_val >= 0.60: return '#FFD700'  # Gold (60-80%)
            else: return '#228B22'  # Green (<60%)

        # Draw transmission elements with loading colors
        if has_lines:
            for line_id, load_val in loading.items():
                line = network.lines.loc[line_id]
                bus0_coords = [bus_coords.loc[line.bus0, 'lat'], bus_coords.loc[line.bus0, 'lon']]
                bus1_coords = [bus_coords.loc[line.bus1, 'lat'], bus_coords.loc[line.bus1, 'lon']]
                color = get_color_by_loading(load_val)
                weight = 0.5 + (load_val * 4)  # Width scales with loading
                folium.PolyLine(
                    [bus0_coords, bus1_coords],
                    color=color,
                    weight=weight,
                    opacity=0.7,
                    tooltip=f'{line_id}: {load_val*100:.1f}% loaded'
                ).add_to(m)
        elif has_internal_links:
            for link_id, load_val in loading.items():
                link = network.links.loc[link_id]
                bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
                bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
                color = get_color_by_loading(load_val)
                weight = 0.5 + (load_val * 4)
                folium.PolyLine(
                    [bus0_coords, bus1_coords],
                    color=color,
                    weight=weight,
                    opacity=0.7,
                    tooltip=f'{link_id}: {load_val*100:.1f}% loaded'
                ).add_to(m)

        # Add small bus markers
        for bus_id, row in bus_coords.iterrows():
            folium.CircleMarker(
                location=[row['lat'], row['lon']],
                radius=2,
                color='orange',
                fill=True,
                fillOpacity=0.6,
                tooltip=bus_id
            ).add_to(m)

        display(m)

        # Print statistics
        print(f'\nTransmission Loading ({component_type}) at {snapshot}:')
        print(f'  Max: {loading.max()*100:.1f}% | Mean: {loading.mean()*100:.1f}%')
        print(f'  >95% loaded: {(loading >= 0.95).sum()} | >80% loaded: {(loading >= 0.80).sum()}')
        print('Color legend: Green (<60%) → Gold (60-80%) → Orange (80-95%) → Red (>95%)')
    else:
        print(f'⚠️  No transmission loading data available at {snapshot}')
except Exception as e:
    print(f'⚠️  Transmission loading map unavailable: {e}')
    import traceback
    traceback.print_exc()
Make this Notebook Trusted to load map: File -> Trust Notebook

Transmission Loading (lines) at 2023-05-18 12:00:00:
  Max: 100.0% | Mean: 8.8%
  >95% loaded: 8 | >80% loaded: 11
Color legend: Green (<60%) → Gold (60-80%) → Orange (80-95%) → Red (>95%)
[7]:
# Interactive network map showing generator dispatch (electricity only)
try:
    import folium

    # Get total generation per bus at snapshot (exclude hydrogen generators)
    h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2']
    elec_gens = network.generators[~network.generators.carrier.isin(h2_carriers)]
    gen_per_bus = network.generators_t.p.loc[snapshot, elec_gens.index].groupby(elec_gens.bus).sum()
    gen_per_bus = gen_per_bus[gen_per_bus > 0]  # Only buses with generation

    if len(gen_per_bus) > 0 and 'x' in getattr(network.buses, 'columns', []):
        # Convert bus coordinates
        t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
        bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
        bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)

        # Create folium map
        center_lat = bus_coords['lat'].mean()
        center_lon = bus_coords['lon'].mean()
        m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')

        # Draw light gray transmission lines in background
        if len(network.lines) > 0:
            for line_id, line in network.lines.iterrows():
                bus0_coords = [bus_coords.loc[line.bus0, 'lat'], bus_coords.loc[line.bus0, 'lon']]
                bus1_coords = [bus_coords.loc[line.bus1, 'lat'], bus_coords.loc[line.bus1, 'lon']]
                folium.PolyLine([bus0_coords, bus1_coords], color='lightgray', weight=0.5, opacity=0.3).add_to(m)

        # Color and size encoding for generation
        max_gen = gen_per_bus.max()
        def color_by_generation(gen_val):
            norm = gen_val / max_gen
            if norm < 0.25: return '#FFE4E1'  # Misty rose
            elif norm < 0.50: return '#FF9999'  # Light red
            elif norm < 0.75: return '#FF6666'  # Red
            else: return '#CC0000'  # Dark red

        # Draw bus markers sized and colored by generation
        for bus_id, gen_mw in gen_per_bus.items():
            coords = bus_coords.loc[bus_id]
            norm_gen = gen_mw / max_gen
            radius = 3 + (norm_gen * 12)  # Size 3-15 pixels
            color = color_by_generation(gen_mw)
            folium.CircleMarker(
                location=[coords['lat'], coords['lon']],
                radius=radius,
                color=color,
                fill=True,
                fillColor=color,
                fillOpacity=0.7,
                tooltip=f'{bus_id}: {gen_mw:.1f} MW'
            ).add_to(m)

        display(m)

        # Statistics
        print(f'\nGenerator Dispatch at {snapshot}:')
        print(f'  Total generation: {gen_per_bus.sum():,.0f} MW')
        print(f'  Active buses: {len(gen_per_bus[gen_per_bus > 0])}')
        print(f'  Max at single bus: {gen_per_bus.max():,.0f} MW')
        print('Color legend: Light gray (no gen) → Light red → Dark red (peak generation)')
    else:
        print('⚠️  No generator dispatch data available at this snapshot')
except Exception as e:
    print(f'⚠️  Generator dispatch map unavailable: {e}')
    import traceback
    traceback.print_exc()
Make this Notebook Trusted to load map: File -> Trust Notebook

Generator Dispatch at 2023-05-18 12:00:00:
  Total generation: 33,422 MW
  Active buses: 2040
  Max at single bus: 3,978 MW
Color legend: Light gray (no gen) → Light red → Dark red (peak generation)
[8]:
# Interactive installed capacity by carrier
generators_p_nom = network.generators.p_nom.groupby(
    network.generators.carrier).sum().sort_values(ascending=True)

# Remove small contributors
generators_p_nom = generators_p_nom[generators_p_nom > 50]

# Create interactive bar chart
fig = go.Figure()

fig.add_trace(go.Bar(
    y=generators_p_nom.index,
    x=generators_p_nom.values,
    orientation='h',
    marker=dict(
        color=[color_map.get(c, '#CCCCCC') for c in generators_p_nom.index],
    ),
    text=[f'{val:,.0f} MW' for val in generators_p_nom.values],
    textposition='auto',
    hovertemplate='<b>%{y}</b><br>Capacity: %{x:,.0f} MW<extra></extra>'
))

fig.update_layout(
    title=dict(text='Generator Capacity by Technology - Historical_2023_etys', x=0.5, xanchor='center'),
    xaxis_title='Installed Capacity (MW)',
    yaxis_title='',
    height=500,
    template='plotly_white',
    showlegend=False
)

fig.show()

Storage Analysis#

[9]:
# Interactive storage dispatch and state of charge visualization
if len(network.storage_units) > 0:
    p_storage = network.storage_units_t.p.sum(axis=1)
    state_of_charge = network.storage_units_t.state_of_charge.sum(axis=1)

    # Create subplot with two y-axes
    fig = make_subplots(
        rows=2, cols=1,
        subplot_titles=('Storage Dispatch', 'State of Charge'),
        vertical_spacing=0.12
    )

    # Storage dispatch
    fig.add_trace(
        go.Scatter(
            x=p_storage.index,
            y=p_storage.values,
            name='Storage Dispatch',
            line=dict(color='#32CD32', width=1.5),
            fill='tozeroy',
            fillcolor='rgba(50, 205, 50, 0.3)',
            hovertemplate='Time: %{x}<br>Power: %{y:.0f} MW<extra></extra>'
        ),
        row=1, col=1
    )

    # State of charge
    fig.add_trace(
        go.Scatter(
            x=state_of_charge.index,
            y=state_of_charge.values,
            name='State of Charge',
            line=dict(color='#00CED1', width=1.5),
            fill='tozeroy',
            fillcolor='rgba(0, 206, 209, 0.3)',
            hovertemplate='Time: %{x}<br>Energy: %{y:.0f} MWh<extra></extra>'
        ),
        row=2, col=1
    )

    fig.update_xaxes(title_text='Time', row=2, col=1)
    fig.update_yaxes(title_text='Power (MW)', row=1, col=1)
    fig.update_yaxes(title_text='Energy (MWh)', row=2, col=1)

    fig.update_layout(
        title=dict(text='Storage Analysis - Historical_2023_etys', x=0.5, xanchor='center'),
        height=700,
        hovermode='x unified',
        template='plotly_white',
        showlegend=True
    )

    fig.show()

    # Storage statistics
    print(f'\nStorage Statistics:')
    print(f'  Total storage capacity: {network.storage_units.p_nom.sum():,.0f} MW')
    print(f'  Total energy capacity: {network.storage_units.max_hours.sum() * network.storage_units.p_nom.sum():,.0f} MWh')
    print(f'  Peak discharge: {p_storage.max():,.0f} MW')
    print(f'  Peak charge: {p_storage.min():,.0f} MW')
    print(f'  Max state of charge: {state_of_charge.max():,.0f} MWh')
else:
    print('No storage units in network')

Storage Statistics:
  Total storage capacity: 5,216 MW
  Total energy capacity: 970,234 MWh
  Peak discharge: 4,366 MW
  Peak charge: -3,978 MW
  Max state of charge: 24,346 MWh

Hydrogen System Analysis#

The hydrogen system models the Power-to-Gas-to-Power pathway:

  • Electrolysis: Converts electricity to hydrogen (stored in central H2 bus)

  • H2 Storage: Central hydrogen storage (copper-plate model)

  • H2 Turbines: Converts hydrogen back to electricity when needed

About the Hydrogen System#

The hydrogen subsystem models the Power-to-Gas-to-Power pathway:

  1. Electrolysis: Converts surplus renewable electricity to hydrogen (dashed turquoise links)

  2. H2 Storage: Stores hydrogen for later use (purple markers)

  3. H2 Turbines: Converts hydrogen back to electricity during high demand (red markers)

This provides seasonal energy storage, shifting summer renewable surplus to winter demand peaks.

[10]:
# Hydrogen System Analysis
# Check if hydrogen system exists
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
has_h2_stores = len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values

if has_electrolysis or has_h2_turbines:
    print('=' * 80)
    print('HYDROGEN SYSTEM SUMMARY')
    print('=' * 80)

    # Electrolysis
    if has_electrolysis:
        electrolysis = network.links[network.links.carrier == 'electrolysis']
        elec_capacity = electrolysis.p_nom.sum()
        elec_consumed = network.links_t.p0[electrolysis.index].sum().sum() / 1000  # GWh
        elec_efficiency = electrolysis.efficiency.mean()
        h2_produced = elec_consumed * elec_efficiency  # GWh H2
        print(f'\nELECTROLYSIS:')
        print(f'  Capacity: {elec_capacity:,.0f} MW ({len(electrolysis)} units)')
        print(f'  Efficiency: {elec_efficiency*100:.0f}%')
        print(f'  Electricity consumed: {elec_consumed:,.1f} GWh')
        print(f'  Hydrogen produced: {h2_produced:,.1f} GWh_H2')

    # H2 Turbines
    if has_h2_turbines:
        h2_turbines = network.links[network.links.carrier == 'H2_turbine']
        turbine_capacity = h2_turbines.p_nom.sum()
        h2_consumed = network.links_t.p0[h2_turbines.index].sum().sum() / 1000  # GWh H2
        turbine_efficiency = h2_turbines.efficiency.mean()
        elec_generated = h2_consumed * turbine_efficiency  # GWh electricity
        print(f'\nH2 TURBINES:')
        print(f'  Capacity: {turbine_capacity:,.0f} MW ({len(h2_turbines)} units)')
        print(f'  Efficiency: {turbine_efficiency*100:.0f}%')
        print(f'  Hydrogen consumed: {h2_consumed:,.1f} GWh_H2')
        print(f'  Electricity generated: {elec_generated:,.1f} GWh')

    # H2 Storage
    if has_h2_stores:
        h2_stores = network.stores[network.stores.carrier == 'H2_gas']
        h2_storage_capacity = h2_stores.e_nom.sum()
        print(f'\nH2 STORAGE:')
        print(f'  Energy capacity: {h2_storage_capacity:,.0f} MWh ({h2_storage_capacity/1000:.0f} GWh)')
        if h2_stores.index[0] in network.stores_t.e.columns:
            soc = network.stores_t.e[h2_stores.index[0]]
            print(f'  State of charge range: {soc.min():,.0f} - {soc.max():,.0f} MWh')
            print(f'  Final state of charge: {soc.iloc[-1]:,.0f} MWh')

    # Round-trip efficiency
    if has_electrolysis and has_h2_turbines:
        roundtrip_eff = elec_efficiency * turbine_efficiency
        print(f'\nSYSTEM METRICS:')
        print(f'  Round-trip efficiency: {roundtrip_eff*100:.1f}%')
        if elec_consumed > 0:
            actual_roundtrip = elec_generated / elec_consumed
            print(f'  Actual energy ratio (elec out / elec in): {actual_roundtrip*100:.1f}%')
else:
    print('No hydrogen system components found in this network.')
    print('(Hydrogen components are added for future scenarios with H2 power generation)')
No hydrogen system components found in this network.
(Hydrogen components are added for future scenarios with H2 power generation)
[11]:
# Hydrogen System Dispatch Visualization
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
has_h2_stores = len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values

if has_electrolysis and has_h2_turbines:
    electrolysis = network.links[network.links.carrier == 'electrolysis']
    h2_turbines = network.links[network.links.carrier == 'H2_turbine']

    # Electrolysis power consumption (positive = consuming electricity)
    elec_power = network.links_t.p0[electrolysis.index].sum(axis=1) / 1000  # GW

    # H2 turbine power generation (p1 is output, negative means generation)
    turbine_power = -network.links_t.p1[h2_turbines.index].sum(axis=1) / 1000  # GW

    # Net hydrogen flow (positive = production, negative = consumption)
    h2_production = elec_power * electrolysis.efficiency.mean()  # GW H2 produced
    h2_consumption = network.links_t.p0[h2_turbines.index].sum(axis=1) / 1000  # GW H2 consumed

    # Create visualization
    fig = make_subplots(
        rows=3, cols=1,
        subplot_titles=(
            'Electrolysis (Electricity Consumed)',
            'H2 Turbines (Electricity Generated)',
            'H2 Storage State of Charge'
        ),
        vertical_spacing=0.1
    )

    # Electrolysis power
    fig.add_trace(
        go.Scatter(
            x=elec_power.index,
            y=elec_power.values,
            name='Electrolysis',
            line=dict(color='#00CED1', width=1.5),
            fill='tozeroy',
            fillcolor='rgba(0, 206, 209, 0.3)',
            hovertemplate='Time: %{x}<br>Power: %{y:.2f} GW<extra></extra>'
        ),
        row=1, col=1
    )

    # H2 turbine power
    fig.add_trace(
        go.Scatter(
            x=turbine_power.index,
            y=turbine_power.values,
            name='H2 Turbines',
            line=dict(color='#FF6347', width=1.5),
            fill='tozeroy',
            fillcolor='rgba(255, 99, 71, 0.3)',
            hovertemplate='Time: %{x}<br>Power: %{y:.2f} GW<extra></extra>'
        ),
        row=2, col=1
    )

    # H2 storage state of charge
    if has_h2_stores:
        h2_stores = network.stores[network.stores.carrier == 'H2_gas']
        if h2_stores.index[0] in network.stores_t.e.columns:
            h2_soc = network.stores_t.e[h2_stores.index[0]] / 1000  # GWh
            fig.add_trace(
                go.Scatter(
                    x=h2_soc.index,
                    y=h2_soc.values,
                    name='H2 Storage',
                    line=dict(color='#9467BD', width=1.5),
                    fill='tozeroy',
                    fillcolor='rgba(148, 103, 189, 0.3)',
                    hovertemplate='Time: %{x}<br>Energy: %{y:.1f} GWh<extra></extra>'
                ),
                row=3, col=1
            )

    fig.update_yaxes(title_text='Power (GW)', row=1, col=1)
    fig.update_yaxes(title_text='Power (GW)', row=2, col=1)
    fig.update_yaxes(title_text='Energy (GWh)', row=3, col=1)
    fig.update_xaxes(title_text='Time', row=3, col=1)

    fig.update_layout(
        title=dict(text='Hydrogen System Dispatch - Historical_2023_etys', x=0.5, xanchor='center'),
        height=900,
        hovermode='x unified',
        template='plotly_white',
        showlegend=True
    )

    fig.show()

    # Correlation analysis
    print('\nDispatch Patterns:')
    print(f'  Peak electrolysis: {elec_power.max():.2f} GW at {elec_power.idxmax()}')
    print(f'  Peak H2 turbine: {turbine_power.max():.2f} GW at {turbine_power.idxmax()}')
    print(f'  Capacity factor (electrolysis): {elec_power.mean() / (electrolysis.p_nom.sum()/1000) * 100:.1f}%')
    print(f'  Capacity factor (H2 turbines): {turbine_power.mean() / (h2_turbines.p_nom.sum()/1000) * 100:.1f}%')
else:
    print('Hydrogen dispatch visualization not available (no hydrogen system components)')
Hydrogen dispatch visualization not available (no hydrogen system components)
[12]:
# Interactive spatial map of hydrogen infrastructure
try:
    import folium

    has_h2_system = (
        (len(network.links) > 0 and any(network.links.carrier.isin(['electrolysis', 'H2_turbine', 'H2 pipeline']))) or
        (len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values)
    )

    if has_h2_system and 'x' in network.buses.columns:
        # Convert all bus coordinates
        t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
        bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
        bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)

        # Create map centered on GB
        center_lat = bus_coords['lat'].mean()
        center_lon = bus_coords['lon'].mean()
        m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')

        # 1. Draw H2 pipelines (H2 pipeline links)
        if len(network.links) > 0:
            h2_pipelines = network.links[network.links.carrier == 'H2 pipeline']
            if len(h2_pipelines) > 0:
                for link_id, link in h2_pipelines.iterrows():
                    if link.bus0 in bus_coords.index and link.bus1 in bus_coords.index:
                        bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
                        bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
                        capacity_mw = link.p_nom
                        folium.PolyLine(
                            [bus0_coords, bus1_coords],
                            color='#00CED1',  # Cyan for H2 pipelines
                            weight=3,
                            opacity=0.7,
                            tooltip=f'H2 Pipeline: {link_id}<br>Capacity: {capacity_mw:.0f} MW'
                        ).add_to(m)
                print(f'Added {len(h2_pipelines)} H2 pipelines to map')

        # 2. Draw electrolysis links (electricity -> H2) and facilities
        if len(network.links) > 0:
            electrolysis = network.links[network.links.carrier == 'electrolysis']
            if len(electrolysis) > 0:
                # Draw electrolysis links connecting electricity bus to H2 storage
                for link_id, link in electrolysis.iterrows():
                    # bus0 is electricity input, bus1 is H2 output/storage
                    if link.bus0 in bus_coords.index and link.bus1 in bus_coords.index:
                        bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
                        bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
                        folium.PolyLine(
                            [bus0_coords, bus1_coords],
                            color='#40E0D0',  # Turquoise for electrolysis links
                            weight=2,
                            opacity=0.6,
                            dash_array='5, 5',  # Dashed line
                            tooltip=f'Electrolysis Link: {link_id}<br>Capacity: {link.p_nom:.0f} MW'
                        ).add_to(m)

                # Draw electrolysis facility markers at electricity input bus
                for link_id, link in electrolysis.iterrows():
                    if link.bus0 in bus_coords.index:
                        coords = bus_coords.loc[link.bus0]
                        folium.CircleMarker(
                            location=[coords['lat'], coords['lon']],
                            radius=6 + (link.p_nom / 500),  # Size by capacity
                            color='#40E0D0',  # Turquoise
                            fill=True,
                            fillColor='#40E0D0',
                            fillOpacity=0.7,
                            tooltip=f'Electrolysis: {link_id}<br>Capacity: {link.p_nom:.0f} MW<br>Electricity bus: {link.bus0}<br>H2 bus: {link.bus1}'
                        ).add_to(m)
                print(f'Added {len(electrolysis)} electrolysis units and links to map')

        # 3. Draw H2 turbine locations (H2 -> electricity)
        if len(network.links) > 0:
            h2_turbines = network.links[network.links.carrier == 'H2_turbine']
            if len(h2_turbines) > 0:
                for link_id, link in h2_turbines.iterrows():
                    # bus1 is electricity output
                    if link.bus1 in bus_coords.index:
                        coords = bus_coords.loc[link.bus1]
                        folium.CircleMarker(
                            location=[coords['lat'], coords['lon']],
                            radius=6 + (link.p_nom / 500),  # Size by capacity
                            color='#FF6347',  # Tomato red
                            fill=True,
                            fillColor='#FF6347',
                            fillOpacity=0.7,
                            tooltip=f'H2 Turbine: {link_id}<br>Capacity: {link.p_nom:.0f} MW<br>Bus: {link.bus1}'
                        ).add_to(m)
                print(f'Added {len(h2_turbines)} H2 turbine units to map')

        # 4. Draw H2 storage locations
        if len(network.stores) > 0:
            h2_stores = network.stores[network.stores.carrier == 'H2_gas']
            if len(h2_stores) > 0:
                for store_id, store in h2_stores.iterrows():
                    if store.bus in bus_coords.index:
                        coords = bus_coords.loc[store.bus]
                        energy_gwh = store.e_nom / 1000
                        folium.CircleMarker(
                            location=[coords['lat'], coords['lon']],
                            radius=8,
                            color='#9467BD',  # Purple
                            fill=True,
                            fillColor='#9467BD',
                            fillOpacity=0.8,
                            tooltip=f'H2 Storage: {store_id}<br>Capacity: {energy_gwh:.1f} GWh<br>Bus: {store.bus}'
                        ).add_to(m)
                print(f'Added {len(h2_stores)} H2 storage facilities to map')

        # Add legend
        legend_html = '''<div style=\"position: fixed; bottom: 50px; right: 50px; width: 240px;
                         background-color: white; border:2px solid grey; z-index:9999; font-size:14px;
                         padding: 10px\">\n
                         <p style=\"margin:0; font-weight:bold;\">Hydrogen Infrastructure</p>\n
                         <p style=\"margin:5px 0;\"><span style=\"color:#00CED1; font-size:20px;\">━━</span> H2 Pipeline</p>\n
                         <p style=\"margin:5px 0;\"><span style=\"color:#40E0D0; font-size:12px;\">- - -</span> Electrolysis Link</p>\n
                         <p style=\"margin:5px 0;\"><span style=\"color:#40E0D0; font-size:20px;\">●</span> Electrolysis</p>\n
                         <p style=\"margin:5px 0;\"><span style=\"color:#FF6347; font-size:20px;\">●</span> H2 Turbine</p>\n
                         <p style=\"margin:5px 0;\"><span style=\"color:#9467BD; font-size:20px;\">●</span> H2 Storage</p>\n
                         </div>'''
        m.get_root().html.add_child(folium.Element(legend_html))

        display(m)
        print('\\n✓ Hydrogen network map displayed')
        print('Marker size indicates capacity (larger = higher capacity)')
    else:
        print('No hydrogen infrastructure found or missing bus coordinates')
except Exception as e:
    print(f'⚠️  Hydrogen network map unavailable: {e}')
    import traceback
    traceback.print_exc()
No hydrogen infrastructure found or missing bus coordinates

Transmission Analysis#

[13]:
# Interactive transmission loading analysis with congestion hotspots
# Handles both lines (ETYS/Reduced) and links (Zonal)
snapshot_idx = len(network.snapshots) // 2
snapshot = network.snapshots[snapshot_idx]

# Determine which transmission elements to analyze
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = len(network.links) > 0

if has_lines:
    # ETYS/Reduced network - use lines
    loading = (network.lines_t.p0.loc[snapshot].abs() / network.lines.s_nom).fillna(0)
    components = network.lines
    element_type = 'Line'
elif has_internal_links:
    # Zonal network - use internal links only (exclude external and hydrogen links)
    h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas']
    is_external = (
        network.links.bus0.str.lower().str.contains('external', na=False) |
        network.links.bus1.str.lower().str.contains('external', na=False)
    )
    is_h2_link = network.links.carrier.isin(h2_carriers)
    internal_links = network.links[~is_external & ~is_h2_link]
    if len(internal_links) > 0 and snapshot in network.links_t.p0.index:
        link_flows = network.links_t.p0.loc[snapshot, internal_links.index].abs()
        loading = (link_flows / internal_links.p_nom).fillna(0)
        components = internal_links
        element_type = 'Link'
    else:
        loading = pd.Series(dtype=float)
        element_type = 'None'
else:
    loading = pd.Series(dtype=float)
    element_type = 'None'

if len(loading) > 0:
    loading_sorted = loading.sort_values(ascending=False)
    top_n = min(20, len(loading_sorted))
    top_elements = loading_sorted.head(top_n)

    # Create labels with bus names
    element_labels = []
    for elem_id in top_elements.index:
        bus0 = components.loc[elem_id, 'bus0']
        bus1 = components.loc[elem_id, 'bus1']
        element_labels.append(f'{bus0} - {bus1}')

    # Color code by congestion level
    colors_list = [
        '#DC143C' if x >= 0.95 else '#FF8C00' if x >= 0.8 else '#FFD700' if x >= 0.6 else '#228B22'
        for x in top_elements.values
    ]

    fig = go.Figure()
    fig.add_trace(go.Bar(
        y=element_labels,
        x=top_elements.values * 100,
        orientation='h',
        marker=dict(color=colors_list),
        text=[f'{val*100:.1f}%' for val in top_elements.values],
        textposition='auto',
        hovertemplate='<b>%{y}</b><br>Loading: %{x:.1f}%<extra></extra>'
    ))

    # Add reference lines
    fig.add_vline(x=95, line_dash='dash', line_color='red', annotation_text='95% (Congestion)')
    fig.add_vline(x=100, line_dash='solid', line_color='darkred', annotation_text='100% (Limit)')

    fig.update_layout(
        title=dict(text=f'Top {top_n} Most Loaded {element_type}s - {scenario_id}', x=0.5, xanchor='center'),
        xaxis_title=f'{element_type} Loading (%)',
        yaxis_title='',
        height=max(400, top_n * 25),
        template='plotly_white',
        showlegend=False
    )
    fig.show()

    # Statistics
    print(f'\n{element_type} Loading Statistics (snapshot: {snapshot})')
    print(f'  Total {element_type.lower()}s: {len(loading)}')
    print(f'  Max loading: {loading.max()*100:.1f}%')
    print(f'  Mean loading: {loading.mean()*100:.1f}%')
    print(f'  {element_type}s >95% loaded: {(loading >= 0.95).sum()}')
    print(f'  {element_type}s >80% loaded: {(loading >= 0.80).sum()}')
else:
    print('No transmission elements found for loading analysis')

Line Loading Statistics (snapshot: 2023-05-18 12:00:00)
  Total lines: 1592
  Max loading: 100.0%
  Mean loading: 8.8%
  Lines >95% loaded: 8
  Lines >80% loaded: 11

Renewable Curtailment Analysis#

[14]:
# Interactive renewable curtailment analysis
renewable_carriers = ['wind_onshore', 'wind_offshore', 'solar_pv']

# Create subplots for each renewable carrier
fig = make_subplots(
    rows=len(renewable_carriers), cols=1,
    subplot_titles=[carrier.replace('_', ' ').title() for carrier in renewable_carriers],
    vertical_spacing=0.08
)

curtailment_stats = {}

for idx, carrier in enumerate(renewable_carriers, 1):
    # Get generators of this carrier
    gens = network.generators[network.generators.carrier == carrier]

    if len(gens) == 0:
        continue

    # Calculate available and dispatched
    p_nom_sum = gens.p_nom.sum()

    # Check if any of these generators have p_max_pu time series data
    # p_max_pu columns are indexed by generator name, not carrier
    gens_with_pmax = [g for g in gens.index if g in network.generators_t.p_max_pu.columns]

    if len(gens_with_pmax) > 0:
        # Calculate available capacity from p_max_pu (weather-dependent availability)
        available = network.generators_t.p_max_pu[gens_with_pmax].multiply(
            gens.loc[gens_with_pmax, 'p_nom'], axis=1).sum(axis=1)
        # Add any generators without p_max_pu at their full capacity
        gens_without_pmax = [g for g in gens.index if g not in network.generators_t.p_max_pu.columns]
        if len(gens_without_pmax) > 0:
            available = available + gens.loc[gens_without_pmax, 'p_nom'].sum()
    else:
        # No time-varying availability - use installed capacity
        available = pd.Series(p_nom_sum, index=network.snapshots)

    dispatched = network.generators_t.p[gens.index].sum(axis=1)
    curtailed = available - dispatched

    # Add available capacity trace
    fig.add_trace(
        go.Scatter(
            x=available.index,
            y=available / 1000,
            name='Available',
            line=dict(color='orange', width=1.5),
            legendgroup=carrier,
            showlegend=(idx == 1),
            hovertemplate='Available: %{y:.2f} GW<extra></extra>'
        ),
        row=idx, col=1
    )

    # Add available capacity trace (orange line showing weather-dependent potential)
    fig.add_trace(
        go.Scatter(
            x=available.index,
            y=available / 1000,
            name='Available',
            line=dict(color='#FF8C00', width=1.5),
            mode='lines',
            legendgroup=carrier,
            showlegend=(idx == 1),
            hovertemplate='Available: %{y:.2f} GW<extra></extra>'
        ),
        row=idx, col=1
    )

    # Add dispatched trace (filled area from zero)
    fig.add_trace(
        go.Scatter(
            x=dispatched.index,
            y=dispatched / 1000,
            name='Dispatched',
            line=dict(color='#228B22', width=1.5),
            fill='tozeroy',
            fillcolor='rgba(34, 139, 34, 0.6)',
            legendgroup=carrier,
            showlegend=(idx == 1),
            hovertemplate='Dispatched: %{y:.2f} GW<extra></extra>'
        ),
        row=idx, col=1
    )

    # Add curtailed as separate trace showing actual curtailment
    curtailed_positive = curtailed.clip(lower=0)  # Only show positive curtailment
    fig.add_trace(
        go.Scatter(
            x=curtailed_positive.index,
            y=curtailed_positive / 1000,
            name='Curtailed',
            line=dict(color='#DC143C', width=1.5, dash='dot'),
            mode='lines',
            legendgroup=carrier,
            showlegend=(idx == 1),
            hovertemplate='Curtailed: %{y:.2f} GW<extra></extra>'
        ),
        row=idx, col=1
    )

    # Update y-axis label
    fig.update_yaxes(title_text='Power (GW)', row=idx, col=1)

    # Calculate statistics
    curtailment_pct = (curtailed.sum() / available.sum() * 100) if available.sum() > 0 else 0
    capacity_factor = (dispatched.sum() / (p_nom_sum * len(network.snapshots)) * 100)

    curtailment_stats[carrier] = {
        'capacity_mw': p_nom_sum,
        'available_mwh': available.sum(),
        'dispatched_mwh': dispatched.sum(),
        'curtailed_mwh': curtailed.sum(),
        'curtailment_pct': curtailment_pct,
        'capacity_factor': capacity_factor
    }

# Update layout
fig.update_xaxes(title_text='Time', row=len(renewable_carriers), col=1)
fig.update_layout(
    title=dict(text='Renewable Curtailment Analysis - Historical_2023_etys', x=0.5, xanchor='center'),
    height=900,
    hovermode='x unified',
    template='plotly_white'
)

fig.show()

# Print statistics
print('\nRenewable Curtailment Statistics:')
print('=' * 80)
for carrier, stats in curtailment_stats.items():
    print(f"\n{carrier.replace('_', ' ').title()}:")
    print(f"  Installed Capacity: {stats['capacity_mw']:,.0f} MW")
    print(f"  Available Energy: {stats['available_mwh']:,.0f} MWh")
    print(f"  Dispatched Energy: {stats['dispatched_mwh']:,.0f} MWh")
    print(f"  Curtailed Energy: {stats['curtailed_mwh']:,.0f} MWh ({stats['curtailment_pct']:.1f}%)")
    print(f"  Capacity Factor: {stats['capacity_factor']:.1f}%")

Renewable Curtailment Statistics:
================================================================================

Wind Onshore:
  Installed Capacity: 12,892 MW
  Available Energy: 970,593 MWh
  Dispatched Energy: 697,733 MWh
  Curtailed Energy: 272,860 MWh (28.1%)
  Capacity Factor: 32.2%

Wind Offshore:
  Installed Capacity: 14,679 MW
  Available Energy: 1,294,766 MWh
  Dispatched Energy: 1,150,376 MWh
  Curtailed Energy: 144,391 MWh (11.2%)
  Capacity Factor: 46.6%

Solar Pv:
  Installed Capacity: 9,019 MW
  Available Energy: 280,556 MWh
  Dispatched Energy: 276,071 MWh
  Curtailed Energy: 4,485 MWh (1.6%)
  Capacity Factor: 18.2%

Advanced Analysis#

Temporal Patterns and System Dynamics#

[15]:
# Interactive heatmap of hourly generation patterns
# Aggregate generation by hour of day and day of year
if len(network.snapshots) > 24:
    total_gen = network.generators_t.p.sum(axis=1)
    total_demand = network.loads_t.p_set.sum(axis=1)

    # Create DataFrame with datetime index
    df_temporal = pd.DataFrame({
        'generation': total_gen,
        'demand': total_demand,
        'hour': total_gen.index.hour,
        'day': total_gen.index.dayofyear
    })

    # Create hourly average pattern
    hourly_pattern = df_temporal.groupby('hour')[['generation', 'demand']].mean()

    # Interactive line plot of hourly patterns
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=hourly_pattern.index,
        y=hourly_pattern['demand'] / 1000,
        name='Average Demand',
        line=dict(color='blue', width=2),
        hovertemplate='Hour: %{x}<br>Demand: %{y:.1f} GW<extra></extra>'
    ))

    fig.add_trace(go.Scatter(
        x=hourly_pattern.index,
        y=hourly_pattern['generation'] / 1000,
        name='Average Generation',
        line=dict(color='red', width=2),
        hovertemplate='Hour: %{x}<br>Generation: %{y:.1f} GW<extra></extra>'
    ))

    fig.update_layout(
        title='Average Hourly Pattern',
        xaxis_title='Hour of Day',
        yaxis_title='Power (GW)',
        template='plotly_white',
        hovermode='x unified',
        height=400
    )

    fig.show()

    print(f'Peak demand hour: {hourly_pattern["demand"].idxmax()}:00')
    print(f'Minimum demand hour: {hourly_pattern["demand"].idxmin()}:00')
else:
    print('Insufficient timesteps for temporal pattern analysis')
Peak demand hour: 11:00
Minimum demand hour: 3:00

Congestion Hotspots Over Time#

Identify when and where transmission congestion occurs.

[16]:
# Calculate transmission loading over time (handles both lines and links)
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = len(network.links) > 0

if has_lines:
    # ETYS/Reduced network - use lines
    loading_pct = (network.lines_t.p0.abs() / network.lines.s_nom * 100).fillna(0)
    components = network.lines
    element_type = 'Line'
elif has_internal_links:
    # Zonal network - use internal links only (exclude external and hydrogen links)
    h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas']
    is_external = (
        network.links.bus0.str.lower().str.contains('external', na=False) |
        network.links.bus1.str.lower().str.contains('external', na=False)
    )
    is_h2_link = network.links.carrier.isin(h2_carriers)
    internal_links = network.links[~is_external & ~is_h2_link]
    if len(internal_links) > 0:
        link_flows = network.links_t.p0[internal_links.index].abs()
        loading_pct = (link_flows / internal_links.p_nom * 100).fillna(0)
        components = internal_links
        element_type = 'Link'
    else:
        loading_pct = pd.DataFrame()
        element_type = 'None'
else:
    loading_pct = pd.DataFrame()
    element_type = 'None'

if len(loading_pct.columns) > 0:
    # Count congested elements at each timestep (>95% loading)
    congested_count = (loading_pct >= 95).sum(axis=1)

    # Create time series plot
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=congested_count.index,
        y=congested_count.values,
        name=f'Congested {element_type}s',
        line=dict(color='#DC143C', width=1.5),
        fill='tozeroy',
        fillcolor='rgba(220, 20, 60, 0.2)',
        hovertemplate='Time: %{x}<br>Congested: %{y}<extra></extra>'
    ))

    fig.update_layout(
        title=f'Transmission Congestion Over Time ({element_type}s >95% Loaded)',
        xaxis_title='Time',
        yaxis_title=f'Number of Congested {element_type}s',
        template='plotly_white',
        height=400
    )

    fig.show()

    # Find most frequently congested elements
    congestion_frequency = (loading_pct >= 95).sum(axis=0).sort_values(ascending=False)
    top_congested = congestion_frequency.head(10)

    if len(top_congested) > 0 and top_congested.iloc[0] > 0:
        print(f'\nMost Frequently Congested {element_type}s:')
        print(f'{element_type}'.ljust(50) + ' Congested Timesteps')
        print('=' * 75)
        for elem_id, count in top_congested.items():
            if count > 0:
                bus0 = components.loc[elem_id, 'bus0']
                bus1 = components.loc[elem_id, 'bus1']
                pct = count / len(network.snapshots) * 100
                print(f'{bus0} - {bus1:<40} {count:>4} ({pct:.1f}%)')
    else:
        print(f'\n✓ No significant transmission congestion detected')
else:
    print('No transmission elements found for congestion analysis')

Most Frequently Congested Lines:
Line                                               Congested Timesteps
===========================================================================
HAMB4A - STAH4A                                    117 (69.6%)
HARK13 - JUNV1-                                    106 (63.1%)
BOLN41 - NINF41                                     85 (50.6%)
WALP41 - RACO41                                     75 (44.6%)
FERO1- - FERO1S                                     71 (42.3%)
FERO1S - INVE1J                                     64 (38.1%)
BHLA11 - GLEN1Q                                     38 (22.6%)
GRNA1- - JUNV1-                                     14 (8.3%)
GLLU1Q - NETS1-                                     12 (7.1%)
INWI1Q - TORN1-                                      7 (4.2%)

Key Performance Indicators#

[17]:
# Calculate KPIs
print(f'SCENARIO: {scenario_id}')
print('=' * 80)
print()

total_demand = network.loads_t.p_set.sum().sum()
total_generation = network.generators_t.p.sum().sum()

# Load shedding
ls_gens = network.generators[network.generators.carrier == 'load_shedding']
load_shedding = network.generators_t.p[ls_gens.index].sum().sum() if len(ls_gens) > 0 else 0

print(f'DEMAND:')
print(f'  Total demand: {total_demand:,.0f} MWh')
print(f'  Load shedding: {load_shedding:,.0f} MWh ({load_shedding/total_demand*100:.2f}%)')
print(f'  Demand satisfied: {total_demand - load_shedding:,.0f} MWh ({(total_demand - load_shedding)/total_demand*100:.2f}%)')
print()

# Renewable share
renewable_gens = network.generators[network.generators.carrier.isin(['wind_onshore', 'wind_offshore', 'solar_pv', 'large_hydro', 'small_hydro'])]
renewable_gen = network.generators_t.p[renewable_gens.index].sum().sum()
renewable_share = renewable_gen / (total_generation + load_shedding) * 100

print(f'GENERATION:')
print(f'  Total generation: {total_generation:,.0f} MWh')
print(f'  Renewable generation: {renewable_gen:,.0f} MWh')
print(f'  Renewable share: {renewable_share:.1f}%')
print()

# System cost
if hasattr(network, 'objective'):
    print(f'OPTIMIZATION:')
    print(f'  Total system cost: £{network.objective:,.2f}')
    if total_generation + load_shedding > 0:
        cost_per_mwh = network.objective / (total_generation + load_shedding)
        print(f'  Cost per MWh: £{cost_per_mwh:.2f}/MWh')
    print()

# Line statistics
max_loading = abs(network.lines_t.p0 / network.lines.s_nom).max().max()
mean_loading = abs(network.lines_t.p0 / network.lines.s_nom).mean().mean()
n_congested = (abs(network.lines_t.p0 / network.lines.s_nom) >= 0.95).sum().sum()

print(f'TRANSMISSION:')
print(f'  Total lines: {len(network.lines)}')
print(f'  Max line loading: {max_loading*100:.1f}%')
print(f'  Mean line loading: {mean_loading*100:.1f}%')
print(f'  Timesteps with congestion (≥95%): {n_congested}')
print()

# Hydrogen system
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values

if has_electrolysis or has_h2_turbines:
    print(f'HYDROGEN SYSTEM:')
    if has_electrolysis:
        electrolysis = network.links[network.links.carrier == 'electrolysis']
        elec_consumed = network.links_t.p0[electrolysis.index].sum().sum()
        print(f'  Electrolysis capacity: {electrolysis.p_nom.sum():,.0f} MW')
        print(f'  Electricity consumed: {elec_consumed:,.0f} MWh')
    if has_h2_turbines:
        h2_turbines = network.links[network.links.carrier == 'H2_turbine']
        elec_generated = -network.links_t.p1[h2_turbines.index].sum().sum()
        print(f'  H2 turbine capacity: {h2_turbines.p_nom.sum():,.0f} MW')
        print(f'  Electricity generated: {elec_generated:,.0f} MWh')
    if has_electrolysis and has_h2_turbines:
        roundtrip = elec_generated / elec_consumed * 100 if elec_consumed > 0 else 0
        print(f'  Round-trip efficiency: {roundtrip:.1f}%')
    print()

print('=' * 80)
SCENARIO: Historical_2023_etys
================================================================================

DEMAND:
  Total demand: 4,805,482 MWh
  Load shedding: 7 MWh (0.00%)
  Demand satisfied: 4,805,475 MWh (100.00%)

GENERATION:
  Total generation: 4,822,551 MWh
  Renewable generation: 2,152,915 MWh
  Renewable share: 44.6%

OPTIMIZATION:
  Total system cost: £141,445,422.73
  Cost per MWh: £29.33/MWh

TRANSMISSION:
  Total lines: 1592
  Max line loading: 100.0%
  Mean line loading: 8.3%
  Timesteps with congestion (≥95%): 631

================================================================================